| Multivariate polynomials in the Horner basis

> {-# LANGUAGE ViewPatterns, BangPatterns #-}

> module Ply (Ply, var, diff, Ply.const,
>             Ply.len, rec_integ, isConst, horner,
>             Named(Named)) where

> import qualified Aux
> import Aux ((!:))
> import qualified Tree
> import Tree (Tree(..))
> import Text.ParserCombinators.Parsec ((<?>),(<|>))
> import Iv (Iv, iv)
> import qualified Iv
> import qualified Text.ParserCombinators.Parsec as Psec
> import qualified Text.ParserCombinators.Parsec.Language as PsecLang
> import qualified Text.ParserCombinators.Parsec.Expr as PsecExpr
> import qualified Text.ParserCombinators.Parsec.Token as PsecTok

> newtype Ply a = P (Tree a) deriving (Eq, Show)

> instance Num a => Num (Ply a) where
>   P (Lf x)       + P (Lf y)       = P $ Lf $! x+y
>   P x            + P (Br [])      = P $ x
>   P (Br [])      + P y            = P $ y
>   P (Br x)       + P (Br y)       = P $ Br $! Aux.zipKeep (|+|) x y
>   P x            + P (Br (y0:y')) = P $ Br (x|+|y0 !: y')
>   P (Br (x0:x')) + P y            = P $ Br (x0|+|y !: x')

>   P (Lf x)       * P (Lf y) = P $ Lf $! x*y
>   P (Lf x)       * P (Br y) = P $ Br $! Aux.map' ((Lf $! x) |*|) y
>   P (Br x)       * P (Lf y) = P $ Br $! Aux.map' (|*|(Lf $! y)) x
>   P (Br [])      * P (Br y) = P $ Br $! []
>   P (Br (x0:x')) * P (Br y) = P $ Br $ Aux.zipKeep (|+|)
>                                              (Aux.map' (x0|*|) y)
>                                              ((Br []) !: (x' .|*|. y))

>   negate (P a) = P $ fmap negate a

>   fromInteger 0 = P $ Br []
>   fromInteger i = P $ Lf (fromInteger i)

>   abs = error "Ply/abs"
>   signum = error "Ply/signum"

> instance Fractional a => Fractional (Ply a) where
>   recip (P (Lf a))   = P $ Lf $ recip a
>   recip (P (Br [a])) = recip (P a)
>   recip _            = error "Ply.recip: division by non-constant"
>   fromRational r     = P $ Lf (fromRational r)

> instance Functor Ply where
>   fmap f (P p) = P $ fmap f p

> horner :: Num a => Ply a -> [a] -> a
> horner (P (Lf c)) _ = c
> horner (P (Br [])) _ = 0
> horner (P (Br (p0:p'))) (x0:x') = horner (P p0) x' + x0 * horner (P (Br p')) (x0:x')
> horner (P (Br (p0:p'))) [] = error "Ply.horner: out of variables"

> const :: a -> Ply a
> const = P . Lf

> (|+|) :: Num a => Tree a -> Tree a -> Tree a
> p |+| q = let P r = P p + P q in r

> (|*|) :: Num a => Tree a -> Tree a -> Tree a
> p |*| q = let P r = P p * P q in r

> (.|*|.) :: (Num a) => [Tree a] -> [Tree a] -> [Tree a]
> p .|*|. q = let P (Br r) = P (Br p) * P (Br q) in r

> var' :: Fractional a => Int -> Iv a -> Tree a
> var' 0     x = Br [Lf $! Iv.lo x, Lf $! Iv.width x]
> var' (i+1) x = Br [var' i x]

> var :: (Fractional a) => Int -> Iv a -> Ply a
> var i x = P (var' i x)


* Differentiation and Differential Recursion

> diffFrst :: Num t => Int -> [Tree t] -> [Tree t]
> diffFrst (i+1) p     = Aux.map' (diffTree i) p
> diffFrst 0     (_:p) = zipWith (\i -> ((Lf $! (fromInteger i)) |*|)) [1 ..] p

> diffTree :: Num t => Int -> Tree t -> Tree t
> diffTree i (Br p) = Br $! (diffFrst i p)
> diffTree i (Lf _) = Br []

> diff :: Num t => Int -> Ply t -> Ply t
> diff i (P p) = P (diffTree i p)

> rec_integ :: Num t => a -> (Ply t -> Ply t -> a -> a) -> Ply t -> a
> rec_integ z _ (P (Br [])) = z
> rec_integ z f (P p) = let d = diffTree 0 p in f (P $ head $ Tree.forest p) (P d) (rec_integ z f (P d))

> len :: Ply t -> Int
> len (P (Lf _)) = 1
> len (P (Br b)) = length b

> isConst :: Ply t -> Maybe t
> isConst (P (Lf c)) = Just c
> isConst _ = Nothing


* Parsing of polynomials with variable names

> newtype Named a = Named ([(String, Iv a)], [Ply a]) deriving (Show)

> instance Fractional a => Read (Named a) where
>     readsPrec _ s = either (\err -> error (show err)) (\x -> [(x,"")]) $ Psec.parse prob "" s 

> l :: PsecTok.TokenParser ()
> l  = PsecTok.makeTokenParser $ PsecLang.haskellDef { PsecTok.reservedOpNames = [".", "*^-", "*^", "*","/","+","-", "<=", "->", "\\/", "^" ] }

> vardec :: Fractional a => Psec.Parser (String, Iv a)
> vardec = do P (Lf a) <- cnst; PsecTok.reserved l "<="; x <- PsecTok.identifier l; PsecTok.symbol l "<="
>             P (Lf b) <- cnst; PsecTok.reserved l "->"
>             return (x, iv a b) <?> "variable declaration"
>          where cnst = ply (Left "No variables allowed here.")

> prob :: Fractional a => Psec.Parser (Named a)
> prob = do PsecTok.whiteSpace l; vars <- Psec.option [] $ Psec.many1 (Psec.try vardec)
>           ps <- Psec.sepBy1 (ply (Right vars)) (PsecTok.reserved l "\\/") ; Psec.eof
>           return $ Named (vars, ps)

> type Vars a = Either String [(String, Iv a)]

> num :: Fractional a => Psec.Parser a
> num = Psec.try (PsecTok.float l >>= return . fromRational . toRational) <|> 
>       (PsecTok.natural l >>= \i -> 
>        Psec.optionMaybe (PsecTok.reserved l ".") >> (return $ fromInteger i))

> brick :: Fractional a => Bool -> Vars a -> Psec.Parser (Ply a)
> brick nums vars = pvar vars <|> 
>                   PsecTok.parens l (ply vars) <|>
>                   (num >>= \c -> if nums then return (P $ Lf c) else fail "Numbers not allowed here.")

> pow :: (Integral a, Fractional a1) => a1 -> a -> a1
> pow a b = if b < 0 then 1 / a^(-b) else a^b

> block :: Fractional a => Bool -> Vars a -> Psec.Parser (Ply a)
> block nums vars = do rcp <- Psec.optionMaybe (PsecTok.reservedOp l "/")
>                      b <- brick (nums || rcp /= Nothing) (maybe vars (\_ -> Left "No variables allowed under '/'.") rcp)
>                      exp <- Psec.optionMaybe (PsecTok.reservedOp l "^" >> PsecTok.natural l)
>                      expTen <- Psec.optionMaybe $ (PsecTok.reservedOp l "*^" >> PsecTok.integer l) <|>
>                                                   (PsecTok.reservedOp l "*^-" >> PsecTok.natural l >>= (return . negate))
>                      return $ maybe id (\_ (P (Lf c)) -> P (Lf (recip c))) rcp 
>                             $ maybe id (flip (^)) exp 
>                             $ maybe id (\e -> (* (P $ Lf $ pow 10 e))) expTen b

> prod :: Fractional a => Vars a -> Psec.Parser (Ply a)
> prod vars = do j1 <- block True vars
>                (block False vars `Psec.chainl1` return (*) >>= return . (* j1)) <|> return j1

> pvar :: Fractional a => Vars a -> Psec.Parser (Ply a)
> pvar (Left m) = fail m
> pvar (Right vars) = do
>   x <- PsecTok.identifier l
>   maybe (fail ("undeclared variable: " ++ x))
>     (\(i, x) -> return $ var i x) (Aux.lookupIx 0 x vars)

> ply :: Fractional a => Vars a -> Psec.Parser (Ply a)
> ply vars = PsecExpr.buildExpressionParser 
>   [ [ PsecExpr.Infix (PsecTok.reservedOp l "*" >> return (*)) PsecExpr.AssocLeft ]
>     , [ PsecExpr.Infix (PsecTok.reservedOp l "+" >> return (+)) PsecExpr.AssocLeft
>       , PsecExpr.Infix (PsecTok.reservedOp l "-" >> return (-)) PsecExpr.AssocLeft
>       , PsecExpr.Prefix (PsecTok.reservedOp l "-" >> return negate) ] ] (prod vars)
